
library(assertthat)
library(reshape2)
library(ggplot2)
## Warning: package 'ggplot2' was built under R version 3.2.3
library(plyr)
parm <- list(h=1e-2, k=10, I=0.6, m=0.001)
##Cap m < 0.2
##cue/h < 0.05
parm$a <- with(parm, 100*h*(k+100)/100)
parm$b <- 5
steady.state <- data.frame(cue=seq(0, 1, length=1e2))
#steady.state$v <- with(parm, b-a*steady.state$cue)
steady.state$v <- with(parm, a*exp(-b*steady.state$cue))
steady.state$check.C <- steady.state$cue*steady.state$v > parm$h
steady.state$check.B <- steady.state$cue*steady.state$v > with(parm, m*k*h/I+h)
summary(steady.state)
## cue v check.C check.B
## Min. :0.00 Min. :0.007412 Mode :logical Mode :logical
## 1st Qu.:0.25 1st Qu.:0.025876 FALSE:9 FALSE:9
## Median :0.50 Median :0.090322 TRUE :91 TRUE :91
## Mean :0.50 Mean :0.221915 NA's :0 NA's :0
## 3rd Qu.:0.75 3rd Qu.:0.315231
## Max. :1.00 Max. :1.100000
assert_that(any(steady.state$check.C & steady.state$check.B))
## [1] TRUE
steady.state$C <- with(parm, h*k/(steady.state$cue*steady.state$v-h))
steady.state$B <- with(parm, steady.state$cue/h*(I-m*steady.state$C))
steady.state$perB <- steady.state$B/(steady.state$C + steady.state$B)
bestCUE <- steady.state$cue[which.max(subset(steady.state, check.C & check.B)$B)]
plotdf <- melt(steady.state, measure.vars=c('C', 'B', 'perB'))
ggplot(subset(plotdf, check.C & check.B)) + geom_line(aes(x=cue, y=value)) + facet_wrap(~variable, scales='free') + geom_vline(xintercept=bestCUE, color='grey')

library(deSolve)
##
## Attaching package: 'deSolve'
##
## The following object is masked from 'package:graphics':
##
## matplot
ans <- data.frame()
parmArr <- data.frame()
for(ii in 1:1e3){
cue <- runif(2)
parm <- list(index=ii, I=0.4, m=0.001,
h1=1e-2, k1=10, cue1=cue[1],
h2=1e-2, k2=10, cue2=cue[2])
parm$a1 <- with(parm, 100*h1*(k1+100)/100)
parm$b1 <- 5
parm$v1 <- with(parm, a1*exp(-b1*cue1))
parm$a2 <- with(parm, 100*h2*(k2+100)/100)
parm$b2 <- 5
parm$v2 <- with(parm, a2*exp(-b2*cue2))
dC <- function(t, y, parms){
B1 <- y[1]; B2 <- y[2]; C <- y[3]
with(parms,{
ans <- c(v1*B1*C/(k1+C)*cue1 - h1*B1,
v2*B2*C/(k2+C)*cue2 - h2*B2,
I-m*C-v1*B1*C/(k1+C)*cue1-v2*B2*C/(k2+C))
return(list(ans))
})
}
timearr <- seq(0, 365*10, length=1000)
evolution <- lsoda(y=c(B1=15, B2=1, C=25), timearr, dC, parms=parm)
temp <- as.data.frame(evolution)
temp$time <- timearr
temp$index <- ii
parm$winner <- names(temp)[2:3][which.max(temp[1000, 2:3])]
if(any(temp[1000, ] < 0)) parm$winner <- NA
parmArr <- rbind.fill(parmArr, as.data.frame(parm))
ans <- rbind.fill(ans, temp)
}
plotdf <- melt(ans, id.vars=c('time', 'index'))
ggplot(plotdf) + geom_line(aes(x=time, y=value, group=index), alpha=0.5) + facet_wrap(~variable, scale='free')

ggplot(parmArr) + geom_point(aes(x=cue1, y=cue2, color=winner)) + geom_abline(slope=1)
